home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
fortran
/
linpkdrv.zip
/
SP.FOR
< prev
next >
Wrap
Text File
|
1984-01-07
|
17KB
|
613 lines
C MAIN PROGRAM
INTEGER LUNIT
C ALLOW 5000 UNDERFLOWS.
C CALL TRAPS(0,0,5001,0,0)
C
C OUTPUT UNIT NUMBER
C
LUNIT = 6
C
CALL SPOTS(LUNIT)
STOP
END
SUBROUTINE SPOTS(LUNIT)
C LUNIT IS THE OUTPUT UNIT NUMBER
C
C TESTS
C SPOCO,SPOFA,SPOSL,SPODI,SPPCO,SPPFA,SPPSL,SPPDI
C SPBCO,SPBFA,SPBSL,SPBDI
C
C LINPACK. THIS VERSION DATED 08/14/78 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C LINPACK SPOCO,SPOSL,SPODI,SPPCO,SPPSL,SPPDI
C LINPACK SPBCO,SPBSL,SPBDI
C EXTERNAL SPOXX,SMACH
C BLAS SAXPY,SDOT,SASUM
C FORTRAN ABS,AMAX1,FLOAT,MAX0
C
C INTERNAL VARIABLES
C
REAL A(15,15),AB(15,15),AINV(15,15),ASAVE(15,15)
REAL AP(120),B(15),SDOT,X(15),XB(15),XEXACT(15)
REAL XP(15),T,Z(15)
REAL ANORM,AINORM,COND,COND1,DET(2),DETB(2),DETP(2)
REAL EN,ENORM,EPS,FNORM,ONEPX,Q(6),QS(6),RCOND,RCONDB
REAL RCONDP,RNORM,S,SASUM,SMACH,XNORM
INTEGER I,INFO,INFOB,INFOP,IQ(6),I1,J,JB
INTEGER K,KASE,KB,KBFAIL,KNPD,KOUNT,KPFAIL
INTEGER KSUSP(6),LDA,LUNIT,M,N,NPRINT
LOGICAL KBF,KPF
C
LDA = 15
C
C WRITE MATRIX AND SOLUTIONS IF N .LE. NPRINT
C
NPRINT = 3
C
WRITE (LUNIT,560)
WRITE (LUNIT,1000)
C
DO 10 I = 1, 6
KSUSP(I) = 0
10 CONTINUE
KNPD = 0
KPFAIL = 0
KBFAIL = 0
C
C SET EPS TO ROUNDING UNIT FOR REAL ARITHMETIC
C
EPS = SMACH(1)
WRITE (LUNIT,570) EPS
WRITE (LUNIT,550)
C
C START MAIN LOOP
C
KASE = 1
20 CONTINUE
C
C GENERATE TEST MATRIX
C
CALL SPOXX(A,LDA,N,KASE,LUNIT)
C
C N = 0 SIGNALS NO MORE TEST MATRICES
C
C ...EXIT
IF (N .LE. 0) GO TO 540
ANORM = 0.0E0
DO 30 J = 1, N
ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1))
30 CONTINUE
WRITE (LUNIT,720) ANORM
C
IF (N .GT. NPRINT) GO TO 50
WRITE (LUNIT,550)
DO 40 I = 1, N
WRITE (LUNIT,760) (A(I,J), J = 1, N)
40 CONTINUE
WRITE (LUNIT,550)
50 CONTINUE
C
C GENERATE EXACT SOLUTION
C
XEXACT(1) = 1.0E0
IF (N .GE. 2) XEXACT(2) = 0.0E0
IF (N .LE. 2) GO TO 70
DO 60 I = 3, N
XEXACT(I) = -XEXACT(I-2)
60 CONTINUE
70 CONTINUE
C
C SAVE MATRIX AND GENERATE R.H.S.
C
DO 90 I = 1, N
B(I) = 0.0E0
DO 80 J = 1, N
ASAVE(I,J) = A(I,J)
B(I) = B(I) + A(I,J)*XEXACT(J)
80 CONTINUE
X(I) = B(I)
XP(I) = X(I)
XB(I) = X(I)
90 CONTINUE
C
C FACTOR AND ESTIMATE CONDITION
C
RCOND = -1.0E0
CALL SPOCO(A,LDA,N,RCOND,Z,INFO)
C
C OUTPUT NULL VECTOR
C
IF (N .GT. NPRINT .OR. INFO .NE. 0) GO TO 110
WRITE (LUNIT,770)
DO 100 I = 1, N
WRITE (LUNIT,780) Z(I)
100 CONTINUE
WRITE (LUNIT,550)
110 CONTINUE
C
C FACTOR PACKED FORM AND COMPARE
C
KPF = .FALSE.
K = 0
DO 130 J = 1, N
DO 120 I = 1, J
K = K + 1
AP(K) = ASAVE(I,J)
120 CONTINUE
130 CONTINUE
RCONDP = -1.0E0
CALL SPPCO(AP,N,RCONDP,Z,INFOP)
IF (INFOP .EQ. INFO) GO TO 140
WRITE (LUNIT,880)
WRITE (LUNIT,920) INFO,INFOP
KPF = .TRUE.
140 CONTINUE
IF (RCONDP .EQ. RCOND) GO TO 150
WRITE (LUNIT,880)
WRITE (LUNIT,930) RCOND,RCONDP
KPF = .TRUE.
150 CONTINUE
K = 0
KOUNT = 0
DO 170 J = 1, N
DO 160 I = 1, J
K = K + 1
IF (AP(K) .NE. A(I,J)) KOUNT = KOUNT + 1
160 CONTINUE
170 CONTINUE
IF (KOUNT .EQ. 0) GO TO 180
WRITE (LUNIT,880)
WRITE (LUNIT,940) KOUNT
KPF = .TRUE.
180 CONTINUE
C
C FACTOR BAND FORM AND COMPARE
C
KBF = .FALSE.
M = 0
DO 200 J = 1, N
DO 190 I = 1, J
IF (ASAVE(I,J) .NE. 0.0E0) M = MAX0(M,J-I)
190 CONTINUE
200 CONTINUE
C
DO 220 J = 1, N
I1 = MAX0(1,J-M)
DO 210 I = I1, J
K = I - J + M + 1
AB(K,J) = ASAVE(I,J)
210 CONTINUE
220 CONTINUE
WRITE (LUNIT,840) M
RCONDB = -1.0E0
CALL SPBCO(AB,LDA,N,M,RCONDB,Z,INFOB)
IF (INFOB .EQ. INFO) GO TO 230
WRITE (LUNIT,830)
WRITE (LUNIT,920) INFO,INFOB
KBF = .TRUE.
230 CONTINUE
IF (RCONDB .EQ. RCOND) GO TO 240
WRITE (LUNIT,830)
WRITE (LUNIT,930) RCOND,RCONDB
KBF = .TRUE.
240 CONTINUE
KOUNT = 0
DO 250 J = 1, N
IF (AB(M+1,J) .NE. A(J,J)) KOUNT = KOUNT + 1
250 CONTINUE
IF (KOUNT .EQ. 0) GO TO 260
WRITE (LUNIT,830)
WRITE (LUNIT,940) KOUNT
KBF = .TRUE.
260 CONTINUE
C
C TEST FOR DEFINITENESS
C
IF (INFO .EQ. 0) GO TO 270
WRITE (LUNIT,860) INFO
KNPD = KNPD + 1
GO TO 530
270 CONTINUE
COND = 1.0E0/RCOND
WRITE (LUNIT,590) COND
ONEPX = 1.0E0 + RCOND
IF (ONEPX .EQ. 1.0E0) WRITE (LUNIT,580)
C
C COMPUTE INVERSE, DETERMINANT AND COND1 = TRUE CONDITION
C
DO 290 J = 1, N
DO 280 I = 1, J
AINV(I,J) = A(I,J)
280 CONTINUE
290 CONTINUE
CALL SPODI(AINV,LDA,N,DET,11)
AINORM = 0.0E0
DO 310 J = 1, N
DO 300 I = J, N
AINV(I,J) = AINV(J,I)
300 CONTINUE
AINORM = AMAX1(AINORM,SASUM(N,AINV(1,J),1))
310 CONTINUE
COND1 = ANORM*AINORM
WRITE (LUNIT,600) COND1
WRITE (LUNIT,800) DET(1)
WRITE (LUNIT,810) DET(2)
C
C SOLVE A*X = B
C
CALL SPOSL(A,LDA,N,X)
C
IF (N .GT. NPRINT) GO TO 330
WRITE (LUNIT,610)
DO 320 I = 1, N
WRITE (LUNIT,790) X(I)
320 CONTINUE
WRITE (LUNIT,550)
330 CONTINUE
C
C MORE PACKED COMPARE
C
CALL SPPSL(AP,N,XP)
KOUNT = 0
DO 340 I = 1, N
IF (XP(I) .NE. X(I)) KOUNT = KOUNT + 1
340 CONTINUE
IF (KOUNT .EQ. 0) GO TO 350
WRITE (LUNIT,880)
WRITE (LUNIT,950) KOUNT
KPF = .TRUE.
350 CONTINUE
CALL SPPDI(AP,N,DETP,11)
IF (DETP(1) .EQ. DET(1) .AND. DETP(2) .EQ. DET(2))
* GO TO 360
WRITE (LUNIT,880)
WRITE (LUNIT,960) DETP
KPF = .TRUE.
360 CONTINUE
KOUNT = 0
K = 0
DO 380 J = 1, N
DO 370 I = 1, J
K = K + 1
IF (AP(K) .NE. AINV(I,J)) KOUNT = KOUNT + 1
370 CONTINUE
380 CONTINUE
IF (KOUNT .EQ. 0) GO TO 390
WRITE (LUNIT,880)
WRITE (LUNIT,970) KOUNT
KPF = .TRUE.
390 CONTINUE
C
C MORE BAND COMPARE
C
CALL SPBSL(AB,LDA,N,M,XB)
KOUNT = 0
DO 400 I = 1, N
IF (XB(I) .NE. X(I)) KOUNT = KOUNT + 1
400 CONTINUE
IF (KOUNT .EQ. 0) GO TO 410
WRITE (LUNIT,830)
WRITE (LUNIT,950) KOUNT
KBF = .TRUE.
410 CONTINUE
CALL SPBDI(AB,LDA,N,M,DETB)
IF (DETB(1) .EQ. DET(1) .AND. DETB(2) .EQ. DET(2))
* GO TO 420
WRITE (LUNIT,830)